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Abstract. We introduce a max-plus analogue of the Petrov-Galerkin finite 
element method to solve finite horizon deterministic optimal control problems. 
The method relies on a max-plus variational formulation. Wc show that the 
error in the sup norm can be bounded from the difl'erence between the value 
function and its projections on max-plus and min-plus semimodules, when 
the max-plus analogue of the stiffness matrix is exactly known. In general, 
the stiffness matrix must be approximated: this requires approximating the 
operation of the Lax-Olcinik semigroup on finite elements. We consider two 
approximations relying on the Hamiltonian. We derive a convergence result, 
in arbitrary dimension, showing that for a class of problems, the error es- 
timate is of order S + Ax((5)~^ or y/S + Ax{&)~^ , depending on the choice 
of the approximation, where <5 and Ax are respectively the time and space 
discretization steps. We compare our method with another max-plus based 
discretization method previously introduced by Fleming and McEneaney. We 
give numerical examples in dimension 1 and 2. 



1. Introduction 
We consider the optimal control problem: 

(la) maximize / ^(x(s), u(s)) ds + 0(x(r)) 

Jo 

over the set of trajectories (x(-),u(-)) satisfying 

(lb) x(.s) = /(x(s), u(,s)), x(s) e X, n{s) G U , 

for all < s < T and 

(Ic) x(0) = X . 

Here, the state space X is a subset of M", the set of control values U is a subset 
of M'", the horizon T > and the initial condition x € X are given, we assume 
that the map u(-) is measurable, and that the map x(-) is absolutely continuous. 
We also assume that the instantaneous reward or Lagrangian £ : X x [/ — > M, and 
the dynamics f : X x U M", are sufficiently regular maps, and that the terminal 
reward </) is a map X ^ M U {— oo}. 
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We are interested in the numerical computation of the value function v which 
associates to any (x,t) £ X x [0,T] the suprcmum v{x,t) of £{x{s) , u{s)) ds + 
(f>{x{t)), under the constraints ljlb|l . for < s < t and ifTcjl . It is known that, under 
certain regularity assumptions, v is solution of the Hamilton- Jacobi equation 

Ov 3v 

(2a) ^—+H{x,—)^Q, (x,t)eXx (0,T] , 

with initial condition: 

(2b) v{x,{)) = (j){x), xeX , 

where H(x,p) = sup^^gj/ £{x, u) + p ■ f{x, u) is the Hamiltonian of the problem (see 
for instance |Lio82llF^lhia,r94j V 

Several techniques have been proposed in the litterature to solve this problem. 
We mention for example finite difference schemes and the method of the vanishing 
viscosity |CL84| , the anti-diffusive schemes for advection |BZ05| , the finite elements 
approach |GR85j (in the case of the stopping time problem) , the so-called discrete 
dynamic programming method or semi-lagrangian method K;D83I, iriDI84| . |Fal87| . 
|FF94j . |FC!99j . |CFF04j . the Markov chain approximations |BD99j . Other schemes 
have been obtained by integration from the essentially nonoscillatory (ENO) schemes 
for the hyperbolic conservation laws (see for instance |OS91p . Recently, max-plus 
methods have been proposed to solve first-order Hamilton- Jacobi equations |MH98j , 

|MH99) . IEHqqI, IHcEoll, IMosl, IM4| . HElol. 

Recall that the max-plus semiring, Mmax, is the set KUj— oo}, equipped with the 
addition a©6 = max(a, b) and the multiplication a®b = a + b. In the sequel, let 5* 
denote the evolution semigroup of (j2J), or Lax-Oleinik semigroup, which associates 
to any map (j) the function := v{-, t), where v is the value function of the optimal 
control problem Maslov |Mas73| observed that the semigroup S"* is max-plus 
linear, meaning that for all maps f,g from X to Mmax, and for all A G Kmax, we 
have 

S*{f(Bg) = S'.f(BS'g , 
S'{\f) = XS'f , 

where f(3g denotes the map x i— > f{x)®g{x), and A/ denotes the map x i— *■ \®f{x). 
Linear operators over max-plus type semirings have been widely studied, see for 
insta nce |rn79llM??92l|Jj(i:iO092|IT?TvT97ll(4Mni^ see also [EtOEl- 

In |FM00j , Fleming and McEneaney introduced a max-plus based discretization 
method to solve a subclass of Hamilton- Jacobi equations (with a Lagrangian (. 
quadratic with respect to u, and a dynamics / afhne with respect to u). They use 
the max-plus linearity of the semigroup S** to approximate the value function by 
a function f ^ of the form: 

(3) vi = sup {A* + Wi} , 

i<i<p 

where {wi}i<i<p is a given family of functions (a max-plus "basis") and {A-}i<i<p 
is a family of scalars (the "coefficients" of on the max-plus "basis"), which must 
be determined. They proposed a discretization scheme in which A* is computed 
inductively by applying a max-plus linear operator to A*~*, where 8 is the time 
discretization step. Thus, their scheme can be interpreted as the dynamic pro- 
gramming equation of a discrete control problem. 
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In this paper, we introduce a max-plus analogue of the finite element method, the 
"MFEM" , to solve the deterministic optimal control problem . We still look for 
an approximation of the form However, to determine the "coefficients" A*, we 
use a max-plus analogue of the notion of variational formulation, which originates 
from the notion of generalized solution of Hamilton-Jacobi equations of Maslov and 
Kolokoltsov KMSai, |KM97I Section 3.2]. We choose a family {zj}i<j<q of test 
functions and define inductively to be the maximal function of the form © 
satisfying 

(4) {vi I zj) < {S'vl' I z,) VI < J < g , 

where (• | •) denotes the max-plus scalar product (see Section |31 for details). We 
show that the corresponding vector of coefficients A* can be obtained by applying to 
A*"'' a nonlinear operator, which can be interpreted as the dynamic programming 
operator of a deterministic zero-sum two players game, with finite action and state 
spaces. The state space of the game corresponds to the set of finite elements. To 
each test function corresponds one possible action of the first player, and to each 
finite element corresponds one possible action of the second player, see Remark |S| 

One interest of the MFEM is to provide, as in the case of the classical finite 
element method, a systematic way to compute error estimates, which can be inter- 
preted geometrically as "projection" errors. In the classical finite element method, 
orthogonal projectors with respect to the energy norm must be used. In the max- 
plus case, projectors on semimodules must be used (note that these projectors 
minimize an additive analogue of Hilbert projective metric |CGQ04| ). 

We shall see that when the value function is nonsmooth, the space of test func- 
tions must be different from the space in which the solution is represented, so 
that our discretization is indeed a max-plus analogue of the Petrov-Galerkin fi- 
nite element method. A convenient choice of finite elements and test functions 
include quadratic functions (also considered by Fleming and McEneaney IFMOOI I 
and norm- like functions, see Sectional 

In the MFEM, we need to compute the value of the max-plus scalar product {z \ 
S^w) for each finite element w and each test function z. In some special cases, {z \ 
S^w) can be computed analytically. In general, we need to approximate this scalar 
product. Here we consider the approximation S^w{x) = w{x) + (5i7(x, Vw(x)), 
for X £ X, which is also used in |MH99| . Our main result, Theorem 1221 provides 
for the resulting discretization of the value function an error estimate of order S + 
Ax{S)~^, where Ax is the "space discretization step", under classical assumptions 
on the control problem and the additionnal assumption that the value function 
is semiconvex for all t £ [0,T]. This is comparable with the order obtained in the 
simplest dicrete dynamic programming method, see |('DI84j . |Fal87| . |('UF89j . To 
avoid solving a difficult (nonconvex) optimization problem, we propose a further 
approximation of the max-plus scalar product (z | S^w), for which we obtain an 
error estimate of order VS + Ax{S)~^, which is yet comparable to the order of the 
existing discretization methods |CDI84j . |Fal87j . ;CDF89j . |CL84j . 

Note that the discretization grid need not be regular: in Theorem Ax is 
defined for an arbitrary grid in term of Voronoi tesselations. 

The paper is organised as follows. In Sectional we recall some basic tools and 
notions: residuation, semimodules and projection. In Section 13 we present the 
formulation of the max-plus finite element method. In Section 01 we compare our 
method with the method proposed by Fleming and McEneaney in [FMOOj . In 
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Sectional we state an error estimate and we give the main convergence theorem. 
Finally, in Sectional we illustrate the method by numerical examples in dimension 
1 and 2. Preliminary results of this paper appeared in |AGL04| . 

Acknov^rledgment: We thank Hcnda El Fekih for advices and suggestions all 
along the development of the present work. 

2. Preliminaries on residuation and projections over semimodules 

In this section we recall some classical residuation results (see for example 
|DJLC53| . |Bir67| . |ET72] . |BCUQ92|), and th eir apphcation to hnear maps on 
idempotent semimodules Csee |LMS01I |CGQ04| ). We also review some results of 
|CGQ96| |CGQ04| concerning projectors over semimodules. Other results on pro- 
jectors over semimodules appeared in |Gon96[ [GM01| . 

2.1. Residuation, semimodules, and linear maps. If {S,<) and {T,<) are 
(partially) ordered sets, wc say that a map f : S T is monotone if s < s' =^ 
/(s) < f{s')- We say that / is residuated if there exists a map f^:T^S such 
that 

(5) m<t^s<fHt). 

The map / is residuated if, and only if, for all t E T, {s E S \ f{s)< t} has a 
maximum element in S. Then, 

f\t) = max{s e S I /(s) <t}, Vt e T . 

Moreover, in that case, we have 

(6) / o /» o / = / and /« o / o . 

In the sequel, we shall consider situations where S (or T) is equipped with an 
idempotent monoid law © [idempotent means that a © a = a). Then the natural 
order on S is defined by a < 6 •^=^ a (B b = b. The supremum law for the natural 
order, which is denoted by V, coincides with © and the infimum law for the natural 
order, when it exists, will be denoted by A. We say that S is complete as a naturally 
ordered set if any subset of S has a least upper bound for the natural order. 

If /C is an idempotent semiring, i.e. a semiring whose addition is idempotent, we 
say that the semiring /C is complete if it is complete as a naturally ordered set, and 
if the left and right multiplications: K, ^ IC, x ^ ax and x i— > xa, are residuated. 
Here and in the sequel, semiring multiplication is denoted by concatenation. 

The max-plus semiring, Kmax, is an idempotent semiring. It is not complete, 
but it can be embedded in the complete idempotent semiring Mmax obtained by 
adjoining +oo to Rmax, with the convention that — oo is absorbing for the multi- 
plication. The map x i— > —x from R to itself yields an isomorphism from Mmax to 
the complete idempotent semiring Mmin, obtained by replacing max by min and by 
exchanging the roles of -l-oo and — oo in the definition of Mmax- 

Semimodules over semirings are defined like modules over rings, mutatis mutan- 
dis, see [l^I SOl , C GQ04| . When /C is a complete idempotent semiring, we say that 
a (right) /C-semimodule X is complete if it is complete as an idempotent monoid, 
and if, for sl\ u € X and A G /C, the right and left multiplications, : X X, 
w I— > wA and : /C ^ A", /i i-^ u/i, are residuated (for the natural order). In a 
complete semimodule X, we define, for all u,v Cz X, 

u\v (L^)Hv) niax{A G IC \ uX < v} . 
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We shall use semimodules of functions: when X is a set and /C is a complete 
idempotent semiring, the set of functions /C^ is a complete /C-semimodule for the 
componentwise addition (u, v) t-^ u(3 v (defined by (u © v){x) = u{x) © ^(2;)), and 
the componentwise multiplication (A, li) 1-^ uX (defined by (uX)(x) = u{x)X). 

If /C is an idempotent semiring, and if X and y are /C-semimodules, we say 
that a map A : X y is linear, or is a linear operator, if for all u,v G X and 
X, ^ G JC, A{uX © vfj.) = A{u)X © A{v)fi. Then, as in classical algebra, we use the 
notation Au instead of A{u). When A is residuated and v £ y, we use the notation 
A\v or A'^v instead of A^{v). We denote by L{X,y) the set of linear operators 
from X to y. If /C is a complete idempotent semiring, if X, y, Z are complete 
A^-semimodules, and if A G L{X,y) is residuated, then L{X,y) and L{X,Z) are 
complete /C-semimodules and the map La '■ L{X,y) L{X,Z), B ^ A o B, is 
residuated and we set A\C := {LAf{C), for aU C G L[X , Z). 

If X and Y arc two sets, /C is a complete idempotent semiring, and a G /C^^^, 
we construct the linear operator A from KX to JC^ which associates to any u G KX 
the function Au G K.-^ such that Au{x) = W y(=Y a{x,y)u{y). We say that A is the 
kernel operator with kernel or matrix a. We shall often use the same notation A 
for the operator and the kernel. As is well known (see for instance |BCOQ92| ), the 
kernel operator A is residuated, and 

iA\v)iy) = A Aix,y)\v{x) . 
In particular, when JC = Mmax, we have 

(7) {A\v){y) = M hA{x, y) + v{x)) = [-A*{-v)]{y) , 

where A* denotes the transposed operator JC^ — > KX , which is associated to the 
kernel A*{ii, x) = A{x, y). (In Q, we use the convention that +00 is absorbing for 
addition.) 

2.2. Projectors on semimodules. Let /C be a complete idempotent semiring and 
V denote a complete subsemimodule of a complete semimodule X, i.e. a subset of 
X that is stable by arbitrary sups and by the action of scalars. We call canonical 
projector on V the map 

(8) Pv.X^X, u 1-^ Pv{u) ^ max{v e V \ V < u} . 

Let W denote a generating family of a complete subsemimodule V, which means 
that any element v gV can be written as w = V{wAu, | w G W}, for some Xnj G /C. 
It is known that 

Pv(u) = V w{w\u) 

(see for instance |CGQ04| ). li B : U ^ X is a residuated linear operator, then 
when U and X are complete semimodules over /C, the image im i? of B is a complete 
subsemimodule of X, and 

(9) Pi„B = BoS« . 

The max-plus finite clement methods relies on the notion of projection on an im- 
age, parallel to a kernel, which was introduced by Cohen, the second author, and 
Quadrat, in |CGQ96| . The following theorem, of which Proposition |21 below is an 
immediate corollary, is a variation on the results of |CG1Q96| Section 6] . 
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Theorem 1 (Projection on an image parallel to a kernel). Let U, X and y he 

complete semimodules over JC. Let B :hi ^ X and C : X ^ y be two residuated 
linear operators over JC Let n§ = B o (C o B)^ o C . We have = o n*^, 
where ^ B o B^ and H'-^ = C" o C. Moreover, is a projector, meaning that 
(ng)^ = ng, and for all x G X : 

ng(x) = max{y e imB \ Cy < Cx} . 

Proof. The first assertion follows from {CoB)^ = B^oCK For the second assertion, 
we have 

(ng)2 = {B o {C o B)^ oC) o {B o {C o bY oC) 
= Bo{CoBfoC (using ®) 

— iig . 

To prove the last assertion, we use that TIb ~ PnaB and (0, we deduce: 

ng(x) = Pi,„BoC«oC(x) 

= max{y £ imB | ?/ < C' o C{x)} 
= niax{?/ G im_B I Cy < Cx} . 

□ 

The results of |CGQ96| characterize the existence and uniqueness, for all x G X, 
oi y E imB such that Cy = Cx. In that case, y = ng(x). 

When /C = Kmax, and C : K-max ^ Il^max is a kernel operator. H*-^ = C" o C has 
an interpretation similar to 

nc^(«) = c« o c{v) = -Pn-^c'i-v) = p-"''^''{v) , 

where — im C* is thought of as a Mmin-subsemimodule of M.^^-^^ and denotes the 
projector on a Mmin-semimodule V, so that, 

p-™'^'(w) = min{u; £ -imC* \ w > v} , 

where < denotes here the usual order on R . When B : M,„„^ Rma^ is also a 
kernel operator, we have 

This factorization will be instrumental in the geometrical interpretation of the finite 
element algorithm. 

Example 2. We take U = {!,■■■ ,p}, X ^ R and Y = {!,■■■ ,q}. Consider the 
linear operators B : ^^^^^^ ^max ^ • ^max ~^ ^max such that 

BX{x) = sup {~-{x ~ Xi f + for aU A e R^a^x , 

l<i<p 2 

and 

{Cf)^ = sup{-a|a; - y^\ + f{x)}, for aU / G if^^ . 

The image of B, imB, is the scmimodule generated in the max-plus sense by the 
functions x t-^ ^^{^ ^ Xif' , for i = 1, • • • ,p. We have 

C'^ IJl{x) = inf {a\x — yi\+ jjii}, for all /i e Rj-,jjj^ , 

l<i<q 
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and the image of C", which coincides with — imC*, is the semimodule generated in 
the min-plus sense by the functions x i-^ a\x — y^j, for i = 1, • • • ,q. 
In fig urc|l(a)| we represent a function v and its projection p-^™'^ {y') (in bold). In 
figure 1(b) I we represent (in bold) the projection Pims(^'~™'^ {v)) — n§(u). 




Figure 1. Example illustrating max-plus and min-plus projectors 



3. The max-plus finite element method 

3.1. M£LX-plus variational formulation. We now describe the max-plus finite 
element method to solve Problem (^. Let V be a complete semimodule of functions 
from X to Mmax- Let 5* : V V and be defined as in the introduction. Using 
the semigroup property 5**+* — S* o , for t, t' > 0, we get: 

(10) v'+^ = S^v\ t = 0,S,--- ,T~6 

with = <j) and 5 = ^, for some positive integer N. Let W C V be a complete 
Mmax-semimodule of functions from X to Kmax such that for all i > 0, u* € W. We 
choose a "dual" semimodule Z of "test functions" from X to Mmax- Recall that the 
max-plus scalar product is defined by 

{u \ v) = sup u{x) + v{x) , 
xex 

for all functions u,v : X ^ Mmax- We replace (|T7)|l by: 

(11) {z\v'+^) ^{z\S^v'), yzez , 

for i = 0, (5, . . . , r - (5, with v\...,v'^ <E W. Equation (jllj can be seen as the 
analogue of a variational or weak formulation. Kolokoltsov and Maslov used this 
formulation in |KM88j and |KM97I Section 3.2] to define a notion of generalized 
solution of Hamilton- Jacobi equations. 
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3.2. Ideal max-plus finite element method. We consider a semimodulc Wh C 
W generated by the family {wi}i<i<p. We call finite elements the functions Wi. 
We approximate v* by S Wh, that is: 

1<1<P 

where A' G Mmax- We also consider a semimodulc Zh <Z Z with generating family 
{zj}i<j<g. The functions 21, • • • , will act as test functions. We replace by 

(12) {z\vl+')^{z\S'vi), ^z^Zn , 

for t = 0, (5, . . . , T — (S, with vf^, . . . ,vj^ £ W^. The function is a given approxi- 
mation of (f). Since Zh is generated hy zi, . . . , Zq, (|12|l is equivalent to 

(13) {z,\vi+')^{z,\S'vi), yi<J<q , 

for i = 0, (5, • • • , T - 5, with u* e W,„ t 0, 5, • • • , T. 

Since Equation p3(l need not have a solution, we look for its maximal subsolu- 
tion, i.e. the maximal solution v*i^^ e Wh of 

(14a) (z, I vi+') < {z, I VI < j < g . 

We also take for the approximate value function dJJ at time the maximal solution 

vl e Wh of 

(14b) vl<v'^ . 

Let us denote by Wh the max-plus linear operator from M^^x ^o W with matrix 
Wh = col(wi)i<i<p, and by the max-plus linear operator from W to Ry^^^ whose 
transposed matrix is Zh = col(zj)i<j<q. This means that WhX ~ \/ i<i<pWiXi for 

all A = {K)i=i,...,p € Kniaxj s-iid {^h'^)j ^ (-^j I ^'o'' all w G yy and j = 1, . . . , (7. 
Applying Theorem ^ to B = Wh and C = Z^ and noting that = im Wh, we 
get: 

Corollary 3. The maximal solution v*^^ G Wh of (|14a|) is given by v*^^ = Sf^{v*f^), 
where 

Note that 11^^ — Pyy^ o p^^h ^ following proposition provides a recursive 
equation verified by the vector of coordinates of u^. 

Proposition 4. Let v\ e Wh be the maximal solution of (|14|l . for t ~ 0, S, . . . ,T . 

Then, for every t — 0,6, ... ,T, there exists a maximal A* £ Kmax such that vj^ ~ 
WhX* , t = 0, 6, ■ ■ ■ ,T, which can be determined recursively from 

(15a) A*+^ - (ZlWh)\{ZlS'WhX') , 

for t — 0, - ■ ■ ,T ~ 6 with the initial condition: 

(15b) A" = Wh\^ . 

Proof. Since G Wh, v]^ = W/jA* for some A* G IRmax and the maximal A* satisfying 
this condition is A* = W^(w^), for all t = 0,6, ... ,T. Since u° is the maximal 
solution of |T4H|) . then by ® and ®, = ^^^,(0) = VF,, o 1^^ (</)), hence A" = 
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Wl o Wh o Wl{<i)) = Wl{(i)). Let i = 5, . . . , T. Using Proposition Theorem □ ® 
and the property that (/ o g)' = o /f for all residuated maps / and g, we get 

A*+* = wlon§^os\Wh\') 

= wl o Wh o wl o {Z*J o Zl o S\WhX') 

- wloiz:)^ozios\WhX') 
= {ziWh)Hzis'WhX*) . 

which yields p5a|l . □ 
For I < i < p and 1 < j < q, we define: 

(16) iMh)J^ - {Zj I w,) 

(17) {Kh)j^ = {zj I S'w,) 

(18) = {iS*)'z, I m,) , 

where S* is the transposed semigroup of S*, which is the evolution semigroup as- 
sociated to the optimal control problem ^ in which the sign of the dynamics is 
changed. The matrices and Kh represent respectively the max-plus linear oper- 
ators Z^Wh and Z^S^Wh- Equation (|15a|l may be written explicitly, for 1 < i < p, 
as 

nrin ( - (A^)^, + max ((A-^)^, + Al.)) • 

1<J<9 ^ l<fe<p / 

Remark 5 . This recursion may be interpreted as the dynamic programming equa- 
tion of a deterministic zero-sum two players game, with finite action and state 
spaces. Here the state space of the game is the finite set {1, • • • ,p} (to each finite 
element corresponds a state of the game). To each test function corresponds one 
possible action j g {1, ■ • ■ ,q} of the first player, and to each finite clement corre- 
sponds one possible action k G {1, • • • ,p} of the second player. Given these actions 
at the state i G {1, • • • ,p}, the cost of the first player, which is the reward of the 
second player, is -~{Mh)ji + {Kh)jk- 

The ideal max-plus finite element method can be summarized as follows: 



Algorithm 1 Ideal max-plus finite element method 

1: Choose the finite elements (ti;i)i<i<p and {zj)i<j<q. Choose the time discretiza- 
tion step S = 

2: Compute the matrix Alh by Hlt)|l and the matrix Kh by H17|) or by (|18|l . 
3: Compute A" = Wh\(f> and = WhX°. 

4: For i = (5, 25, . . . , T, compute A* = Mh\{KhX*-^) and vl = WhXK 



Remark 6. Since w*^ G Wh Vt = 0, • • • , T, the dynamics of w^can be written as a 
function of the matrices Mh and Kh'- 

(19) vi+' = WhoMioKhoWl{vl) . 
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3.3. EfTective max-plus finite element method. In order to implement the 
max-plus finite clement method, we must specify how to compute the entries of the 
matrices A'h and Kh in lfTC|) and ((T7|l or lfTH|) . 

Computing from H16() is an optimization problem, whose objective function is 
concave for natural choices of finite elements and test functions (see Section [S] be- 
low). This problem may be solved by standard optimization algorithms. Evaluating 
every scalar product {z \ S^w) leads to a new optimal control problem since 

(z I S^w) = max 2(x(0)) + / £(x(s), u{s))ds + w(x(5)) , 

Jo 

where the maximum is taken over the set of trajectories (x(-), u(-)) satisfying ljlb|l . 
This problem is simpler to approximate than Problem Q), because the horizon S is 
small, and the functions z and w have a regularizing effect. 

We first discuss the approximation of S^w for every finite element w. The Hamilton- 
Jacobi equation (Pa|) suggests to approximate S^w by the function [S^w]h such that 



(20) [S^w]h{x) ^ w{x) + SH{x, Vw{x)), for all xe X . 

Let [S^Wh]H denote the max-plus linear operator from Mf^ax to W with matrix 
[S^WhjH = col([5"^u'i]/f)i<i<p, which means that 

[S'Wh]H\^ V [S'w,]HX^ 

l<i<'P 

for all A ~ {Xi)i<i<p G K^iax- The above approximation of S^w yields an approx- 
imation of the matrix K^. by the matrix Kjf ji ;= Z^[S'''W^/,]_ff , whose entries are 
given, for 1 < i < p and 1 < j < q, by: 

(21) {KH.h)]i = snp{zj{x) +Wi{x) +6H{x,'Vwi{x))) . 

Thus, computing Khji requires to solve an optimization problem, which is nothing 
but a perturbation of the optimization problem associated to the computation of 
M . We may exploit this observation by replacing Ku.h by the matrix KH,h with 
entries 

(22) {KH,h)ji = {zj\wi)+5 sup H{x,\/wi{x)) , 

for 1 < j < p and 1 < J < Here, argmaxjzj + w,} denotes the set of x such 
that Zj{x) + Wi{x) = {zj I Wi). When this set has only one clement, 1)22(1 yields a 
convenient approximation of Kh- 

Of course, Wi must be differentiable for the approximation H2()|l to make sense. 
When Wi is iion-differentiable, but Zj is differentiable, the dual formula H18|l suggests 
to approximate {Kh)ji by 

sup(zj(a;) + SH{x, —Vzj{x)) + Wi{x)) . 
We may also use the dual formula of (f^ . where S/wi{x) is replaced by — Vzj(x). 
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4. Comparison with the method of Fleming and McEneaney 

Fleming and McEneaney proposed a max-plus based method [FMOOj , which also 
uses a space Wh generated by finite elements, wi, . . . ,Wp, together with the linear 
formulation H10(l . Their method approaches the value function at time t, w*, by 
Wh^J-^ , where Wh = col(wi)i<i<p as above, and is defined inductively by 

(23a) / = 

(23b) = {Wn\{S'WH))^i' , 

for < = 0, (5, . . . , T— 5. This can be compared with the limit case of our finite element 
method, in which the space of test functions Zh is the set of all functions. This 
limit case corresponds to replacing by the identity operator in H15a|) . so that 

(24) A*+^ = Wh\{S^WhX') . 

Proposition 7. Let (/i*) be the sequence of vectors defined by the algorithm of 
Fleming and McEneaney, (|23(l ; let (A*) be the sequence of vectors defined by the 
max-plus finite element method, in the limit case (|24() ; and let w* denote the value 
function at time t. Then, 

(25) Wh^i^ < WhX* <v^ , for < = 0, (5, . . . , T . 

Proof. We first prove that WhX* < f * for i — 0, (5, . . . , T. This can be proved by 
induction. For t = we have WhX° < w° by iMbll . We assume that W^ A* < w*. 
Using (|24|l . we have 

= n^^{s'{WhX*)) . 

Using the monotonicity of the semigroup , we obtain 

WhX'+' < n^^(^^*) 



The second inequality is also proved by induction. For t = 0, we have = A" = 
Wh\<^- Suppose that fi* < A*. By definition of Wh\{S^Wh), we have 

Wh{Wh\S^Wh) < S^Wh , 



hence 



< 



<5w \t 



< S'WhX 



Since 



1*4-5 _ TiA \ ( qS 



Wh\{S'WhX') 



= max{A e RLx I W^X < S'WhX*} , 

we get that fi*+^ < A*+'^. Then yu* < A* for t = 0,6, ...,T. Since Wh is monotone, 
we deduce (Esl. □ 



An approximation of (|23b() using formulae of the same type as H20|l is also dis- 
cussed in |MH99j . 
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5. Error analysis 

5.1. General error estimates. In the sequel we denote by ||u||oo = sup^^j ^ 
M U {+00} the sup-norm of any function v : I ^ W. We also use the same no- 
tation ||w||oo = supjgj \vi\ for a vector v ~ [vi)i£i. For any two sets / and J, 
a map <I> : ^ R"' is said monotone and homogeneous if it is monotone for 
the natural order and if for all u e and A G K, $(u + A) = + A with 

{u + A)(z) = u{i) + A. Monotone homogeneous maps are nonexpansive for the 
sup- norm: — $( see |CT80j . In particular, max-plus 

or mill-plus linear operators are non-expansive for the sup-norm. This property 
will be frequently used in the sequel. In order to simplify notations, we denote 
= {0, 5, • • • , T}, T+ = fs\{0} and r," = fs\{T} . 

Remark 8. To establish the main result of the paper (Theorem |23 below) . we 
shall need only to take the norm of finite valued functions. However, we wish 
to emphasize that all the computations that follow are valid for functions with 
values in R if one replaces every occurence of a term of the form ||u — w|jtx3 by 
doo{u, v) = inf{A > | —X+v < u < X+v}. Observe that doo{u, v) is a semidistance 
and that dooiu, v) ~ \\u — t;||oo, ii u — v takes finite values. Observe also that if a 
map <i> : M R is monotone and homogeneous, doo{^{u),^{v)) < doo{u,v), for 
all u, e M . 

The following lemma shows that the error of the ideal max-plus finite element 
method is controlled by the projection errors ||n^^(u*) — w*||oo- This lemma may 
be thought of as an analogue of Cea's lemma in the classical analysis of the errors of 
the finite element method. Projectors over semimodules in the MFEM correspond 
to orthogonal projectors in the classical finite element method. 

Lemma 9. For t G f^, let u* he the value function at time t, and vj^ be its approx- 
imation given by the ideal max-plus finite element method. We have 



(26) H - ^'"^lloo < lln^Jz;") - ^;"||oo + ^ ~ «*l 



00 



Proof. For all t E Tg , we have 

H+'-v'+'\\oo < \H+'-siiv')u + \\si{v')~v'+'\\^ 

Since is a non-expansive operator, we deduce 

H^' «*+'l|oo < H - «1oo + l|n^^(z;*+^) - «*+^|U . 

The result is obtained by induction on i, using the fact that = PvVdl'f^'^) = 

To obtain an error estimate, we need to bound ||n^^(w*) — w*||oo for all t G . 
Since = o 11^'' , we have 

l|ng.(^*)-^;*lloo = ||^^.^o^^^(^;*)-^,*|u 
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and since ^ non-expansive operator, we get 

(27) ~ «*iu < m'Hv') - v'\u + mwM) - ^''lu • 

Using this inequality together with Lemma |3 we deduce the fohowing corollary. 

Corollary 10. For t € fs, let w* be the value function at time t, and vj^ be its 

approximation given by the ideal max-plus finite element method. We have 

- tets 



hi - v'^Woo < (1 + ^)( sup(||^^':(^;*) - z;*|U + ^wM) - ^;*|U)) 



The following general lemma shows that the error of the effective finite element 
method is controlled by the projection errors and the errors resulting from the 
approximation of the matrix Kh by a matrix Kh. 

Lemma 11. For t €z fg, let be the value function at time t, and v^^ be its 
approximation given by the effective max-plus finite element method, where Kh is 
approximated by Kh- We have 

hi - v-^Woo < (1 + ^)( sup(||n^'^(«*) - «*|U + mwM) - ^'*lloo) 

+ \\kh - KhWoo 

Proof. Since is computed with the approximation Kh of A'^, we have w*^ = W/iA*, 
t G fs, with 

A*+^ = Ml o {KhX') = Wl o {Z*J o {KhX') . 

We have 

ll<+'-^*+'lloo < H^'~sU\\^ + \\sivi-sy\\^ + \\sy-v'+'\\^ 



< 



/l 1 1 OO ) 5 



+hi~v'\\oo + m§,y^')-^'^'\\oo 

< \\KhX' - KhX'Woo + H~ v'\U + m^liv'^') - v'+'\\oo 

< max_^ \{Khh. - {KhU + H " «*lloo + \\K\iv'^') - v'+'\U 

l<i<p 

We deduce that 

hi - v'^Woo < mwjv') - v'\U. + J2 - «*lloo + \\Kh K, 

and so 

T 

)\ 

- tefs 

+ \\Kh-Kh 

□ 

Corollary 12. For t G t^, let w* be the value function at time t, and v\ be its 
approximation given by the effective max-plus finite element method, implemented 



hi - v'^Woo < (1 + ^)( sup(||n^':(«*) - «*iu + mwM) - v^u 
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with the approximation Khm of Kh, given by H21|) . We have 

WvJ: - v^Woo < (1 + ^)( sup(||^^':(^;*) - «*||oo + W^wM) ^ 



max \\[S^w,]h-S' 



i5„ 



l<i<p 



2 1 1 oo 



Proof. Using the same technique as in the precedent lemma and using that Kh^h = 

WKn.h - KhW^ < WiS^WhU - S^WhWoo 
(28) = max \\[S^w,]h ~ S^w,\\oo , 

l<i<p 

which ends the proof. □ 

Corollary 13. For t Cz fs, let w* be the value function at time t, and vj^ be its 
approximation given by the effective max-plus finite element method, implemented 
with the approximation Kh.h of Kh, given by We have 

\\vj: - v-^Woo < (1 + |)f supdin^'^z;* - «*iu + m^-y - v^u 



tefs 



max II [5 Wi]H - S w.^W^o + \\KH,h - KH,h\\oo 

l<i<p 



Proof. We use Lemma together with Equation (|28|l and 

\\KH,h - KhWoo < \\KH,h - KhmWoo + WKhm - Kh\ 



□ 



5.2. Projection errors. In this section, we estimate the projection errors resulting 
from different choices of finite elements. Recall that a function / is c-semiconvex 
if f{x) + infill, where || • II2 is the standard euclidean norm of R", is convex. A 
function / is c-semiconcave if — / is c-semiconvex. Spaces of semiconvex functions 
were intensively used in the max-plus based approximation method of Fleming and 
McEnea n ey |FMnnj , see also |MH98j . |MH99j . |McEn2j . |McEn8| . jMcEn4j . pl87j . 
EdTH, |(TI)K89|. 

We shall use the following finite elements. 

Definition 14 (Pi finite elements). We call Pi finite element or Lipschitz finite 
element centered at point x £ X, with constant a > 0, the function w{x) = 
— a||a; — a;|| 1 where ||x|| 1 = X]"=i 1^*1 ^^e /-'--norm of R". 

The family of Lipschitz finite element of constant a generates, in the max-plus 
sense, the semimodule of Lipschitz continuous functions from X to R of Lipschitz 
constant a with respect to || • ||i. 

Definition 15 (P2 finite elements). We call P2 finite element or quadratic finite 
element centered at point x & X, with Hessian c > 0, the function w{x) = —^\\x — 

When X — R", the family of quadratic finite elements with Hessian c generates, 
in the max-plus sense, the semi-module of lower-scmicontinuous c-semiconvex func- 
tions with values in R. 

Notations. Let F be a subset of R" and / be a function from Y to R. We will 
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denote by Convy the convex hull of Y, riY the relative interior of Y, dom/ the 
effective domain of / and df{x) the subdifferential of / at 2; G dom/. 
When C is a nonempty convex subset of R" and c > 0, a fonction is said to be 
c-strongly convex on C if and only if / — ^c|| • is convex on C . A function / is 
c-strongly concave on C if — / is c-strongly convex on C. 

Let P be a finite subset of M". The Voronoi cell of a point p S P is defined by 

V{p) = {x e M" I \\x- ph < \\x - qhyq e P}. 

The family {V{p)}p^p constitutes a subdivison of M", which is called a Voronoi 
tesselation (see |SUOO| for an introduction to Voronoi tesselations). We define the 

restriction of V{p) to X to be: 

Vx{p) = v{p)nx. 

We define px{P) to be the maximal radius of the restriction to X of the Voronoi 
cells of the points of P: 

px(^') := sup sup ||a;-p||2. 

p6-Px6V0f (p) 

Observe that 

Px{P) sup inf ||a; -p||2. 

The previous definitions are illustrated in Figure [5| The set X is in light gray, 
P ~ {pi, ■ ■ ■ ,pio}, Vx{Pg) is in dark gray and px{P) is indicated by a bidirectional 
arrow. 




Figure 2. Voronoi tesselation 

The next two lemmas bound the projection error in term of the radius of Voronoi 
cells. 

Lemma 16 (Primal projection error). Let X be a compact convex subset of R" . 
Let w : X — > M 6e c-semiconvex and Lipschitz continuous function with Lipschitz 
constant L^ with respect to the euclidean norm. Let Vc(x) = v{x) + |||a;||2. Let 
X ~ ^ + 62(0, ■^), let Xh be a finite subset ofW^ and let Wh denote the complete 
subsemimodule o/R^jj^ generated by the family {wxf^)j.^^x^ where Wx^ix) = ~^\\x— 
XhWl- Then 

\\v - PwiMloo < cdiamXpj^{Xh) 
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Proof. Let W denote the complete subsemimodule of K^^ax generated by the family 
(""^xj^gx- We will first prove that for all x € X, Pwv{x) — v{x). It is obvious that 
Vx e X, v{x) > Pwv{x). Using that Pyv = W o W^, with W = col(w£)^g^' "^^ 
obtain 

Pwv{x) = sup (-^\\x-^l+ M^{^\\y-x\\l+v{y)) ^ 



xex 



sup [ - ^||a;||2 + cx ■ X - sup (cx ■ y - Vc{y)) 



xex ^ 2 yex 

= -■^l|a;||2 + sup (cx • a; - w*(ci)) , 
^ xex 

where v* denotes the Fenchel transform of Vc- Since Vc is l.s.c, convex and proper, 
we have for all a; G X 

(29) Vcix) = v2*{x) = sup {9-x- v*{9)) . 

Using Theorem 23.4 of |R,oc70j . for all x E ri(domwc), the subdifferential of Vc at 
X, dvc{x) = {6 € R"- \ Vc{y) — Vc{x) > 9 ■ {y - x),\/y S X}, is non-empty. Then 
9 G dvc{x) if and only \i v*{9) = 9 ■ x — Vc{x) and consequently, the supremum 
of H29|) is attained for all elements 9 of dvc{x). 

Set q{x) — §\\x\\l. Using the fact that q{y) — q{x) = q'{x) ■ {y - x) + 0{\\y — xWl) 
and that v is Lipschitz continuous with Lipschitz constant L^, we obtain dvc{x) C 
B2{cx, Ly) for all x G liX. Therefore, for all x E liX, 

(30) Vc{x) = sup (ci • X — v*{cx)) . 

xex 

By continuity in the members of Equation H30() . we have the equality for all x G X, 
and so 

Pyvv{x) = —^\\x\\l + SUp(^CX-X — V*{cx)) 

2 xex 

= -^\\x\\l+Vc{x) 

= v{x) , 

for all X € X. 

Now, fix X G X. For x E X, we set ip{x) = cx ■ x — v* {cx). Since P-w^v < Pwv < v, 
we have for all x e X 

Q < v{x) ~ Pwh.v{:x) = Pwv(x) - PwhV{x) 

= sup (p{x) — sup ip{xh) 
xex xheXh 

= sup inf ip{x) — ip{xh) . 
xex SfiSX,, 

We have d{—ip){x) = —cx + cdv*{cx). Since dv* C X, we have d{—ip){x) C c{X — 
x) C 62(0, cdiamX). Hence, ip is Lipschitz continuous with Lipschitz constant 
~ cdiamX. Then for all x E X 

v{x) - Pwh^ix) < sup inf L^\\x~Xh\\2 
xex^h&Xh 

— c diam. X pj^{Xh) ■ 

□ 
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Lemma 17 (Dual projection error). Let X be a bounded subset of M" and X 
a finite subset of M". Let v : X ^ R be a given Lipschitz continuous function 
with Lipschitz constant Ly with respect to the euclidean norm. Let Zj^ denote the 
complete semimodule ofM.^^^ generated by the Pi finite elements {zx)j.^x centered 
at the points of X with constant a > Ly. Then 

\\v - P-(^-*)«||oc < nia + Ly)px{X). 

Proof. It is clear that p-(^x)v > v and using that p-(^x) = [Z*)^ o Z*, with 
Z = col(zx:) , we obtain 

P^'-^x)y(^x) — v{x) ~ inf ( a||a; — a;||i + sup ( — a||y — +«(?/)— ) , 
xex ^ yex ' 

for all X ^ X. Since v is L^j-Lipschitz continuous, we have 



< 


inf 1 


\a\ 


\x — 


x\ 


1 " 


V sup ( 


-a\ 


\y 




SGX 










yex 






< 


inf 1 


{a\ 


\x- 


x\ 


ll ' 


V sup ( 


-a\ 


\y 














yex 






< 


inf 1 


(a\ 


\x — 


x\ 


1 " 


V sup ( 


-a\ 


\y 














yex 







111 + Ly\\y - x\\ 

+Ly\\x - x||i)^ 

= inf ((a + L„)||.T - .t||i + sup(Pu - a)||?/ - .t||i ) 
xex ^ yex ' 

Since a> Ly^ wc deduce 

p-'^^x)v{x) - v{x) < {a + Ly) svip inf ||a; - .xlli < 7i(a + ■ 

xex xex 



□ 



5.3. The approximation errors. To state an error estimate, we make the fol- 
lowing standard assumptions (see |Bar94j for instance): 

- {H\) f : X X U ^ R" is bounded and Lipschitz continuous with respect to 
X, meaning that there exist Lf > and Mf > such that 

\\f{x,u)- f{y,u)\\2 <L/||a;-?;||2 yx,y<EX,u£ U, 

\\f{x,u)\\2<Mf, VxeX^ueU . 

- {H2) ^ : X X J7 — > M is bounded and Lipschitz continuous with respect to 
x, meaning that there exist Lg > and ALg > such that 

\eix, u) - £{y, u)\ < Le\\x - y\\2 Vx, y G X, w G U, 

\£ix,u)\ <Mi, Wx e X,u eU . 

5.3. L Approximation of S^w. 

Lemma 18. Let X be a convex subset o/M". We make assumptions (HI) and (H2). 
Let w : X —* M. be such that w is on a neighborhood of X , Lipschitz continuous 
with Lipschitz constant Ly, with respect to the euclidean norm, ci-semiconvex and 
C2-semiconcave. Then there exists Ki > such that \\[S^w]h — S^w\\od < KiS^, 
for S > 0, where [S^uj]h is given by (|20|l . 
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Proof. We first show that there exists Ki > sucli that 

[S^w]h{x) - S^w{x) > -Ki6^, VxeX . 

For all X Cz X and u Cz U, define :x.u,x to be the trajectory such that x„ 2,(s) = 
.fi^u,xis),u) and x„2,(0) = x. In other words, we apply a constant control u. We 
have 

(S"'w)(a:) > sup{ / £(x„,^(s), M)ds + w(x„,x('5)) | u e U} . . 



Since £ is Lipschitz continuous and / is bounded, we have 



[£{-Xu,x{s),u) - £{x,u)]ds < Li \\:x.u,x{s) - x\\2ds 



< Lf, I Mfsds 



then 
(31) 

Therefore 



[£{;Ku,x{s),u) -£{x,u)]di 



iS^w)ix) > -^LeMfS^ + snp{6£{x,u) + w{xu^x{S)) \ u e U} . 
Since w is Lipschitz continuous and / is bounded and Lipschitz continuous, we have 
w{:x.u,xiS)) - w{x + 6f{x,u)) < Lw\\xu^xiS) ~ x - Sf{x,u)\\2 



< U 



\f{yiu,x{s),u) - f{x,u)\\2ds 



i/l|x„,j;(s) - x\\2ds 



< L^Lf I Mfsds 



'ft 



and so 
(32) 



w{yiu,x{5)) - w{x + Sf{x, u)) 



Moreover, since w is ci-semiconvex, we have 

(33) w{x + Sf{x, u)) > w{x) + SVw{x) ■ f{x, u) - ^M^S^ . 

We deduce from ^ and ^ 

iS^w){x) > -{L(,Mf + L^LfMf+ciM])—+w{x) 

+ sup {S£{x,u) + SWw{x) ■ f{x,u)} 

ueu 

§2 

+ i(;(a:) + (5i7 (x, Vu;(a;)) . 

This ends the first part of the proof. 

We now prove an opposite inequality. For all a; € X and for all measurable 
functions u : [0,(5] — > C/, define Xu,a; to be the trajectory such that Xu,a;(s) = 
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/(xu,x(s), u(s)) and Xu,x(0) = x. Since i{x,u) < H{x,p) —p - f{x,u), for all 
p € R", X £ X and u £ U, we deduce that 

iS^w){x) < sup { ^ iJ(xu,^(s), Vu;(a;))ds + u;(xu,^(J)) 

- \7w{x) ■ /(xu,x(s), u(s))ds I u : [0, (5] ^ [/} 

= sup|^ i^(xu,a;(s), Vu;(a;))ds 

+ w(x(5)) - V?«(a;) • (xu,.((5) - x) | u : [0, (5] ^ L/} . 

Using the fact that t and / are Lipschitz continuous with respect to cc, we have for 
all x,x' £X,p£ M" 

H{x,p)-H{x\p)\ < {Le + Lf\\p\\2)\\x-x'\\2 , 

therefore 

iS^w){x) < sup|(i(? + f \\x^^^{s)^ x\\2ds + 5H{x,Vw{x)) 



+u'(xu,x(<5)) - S/w{x) ■ (xu,,(5) - x) I u : [0,(5] ^ u} 



< {Le + LfLy,)Mf — + 5H{x,yw{x)) 

+ SUp|w(Xu,:r(<5)) - Vw(x) • (Xu,a;((5) - x) | U ! [0, (5] ^ J/j . 

Since w is C2-scniiconcave, we have 

u'(xu..(<5)) < u;(.t) + Vw{x) ■ (xu,,(J) - x) + yM^-J^ ^ 

We obtain 

iS^w){x) < {Li, + LfMD^ + c2Mf)Mf— + w{x) + SH{x,Vw{x)) . 
To end the proof, we take i^i = 5 (LgMj + L fMowM / + max(ci , C2)M|) . □ 
5.3.2. Approximation of the matrix Kh by the matrix Kh- 

Lemma 19. Let X be a compact subset o/R". We consider an upper semicon- 
tinuous function ip : X W and a Lipschitz continuous function tJj : X —i- R with 
Lipschitz constant L^ with respect to a norm || • ||. For e > 0, we define: 

(34a) ^ {x e X \ ip{x) > sup Lp{x') - e}, 

x'ex 

(34b) 5(e) = sup d{x,Fo) , 

xeF, 



where d{x, Fq) ~ infygi?^ \\y — x\\ . We have: 

I sup (<p(x) + 6ip{x)) — [ sup ifi{x 
xgx xex 

where M = supj.g_y ipix) — inf-^gx ipix) 



I sup (<p(x) + (5'0(x)) — [ sup (p{x) + S sup 'fpix)] \< L,p5g[5M) , 

xGX xeX xSargmaxvJ 
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Proof. Since is u.s.c. and X is compact, Fq ~ argmaxi^ and 
(35) sup {f(x) + 6ip{x)) > sup tfiix) + 6 sup ip{x) . 

For e > we have: 

sup (1^9(2;) + Siplx)) = max [ sup (fix) + S^p{x)), sup (fix) + Stp{x))] . 
xex xeF^ xex\F^ 

Let e = (5(sup2.gjf 'ipix) — iiiix^x ^j{x)) = MS (which is finite since ip is continuous 
and X is compact). We have: 

sup {fix) + Sip{x)) < —e+ sup ip{x) + 6 sup ip{x) 

x£X\F^ x£X x£X 

= sup (p{x) + S inf 'il'{x) 
xGF^ ^ex 



< 



sup [(p{x) + Sip{x)'\ 
xeFe 



Therefore 



sup {}pix) + S^{x)) = sup {ip{x) + S4'{x)) 

xGX xS-Fe 



(36) < sup (p{x) + S sup ipi^) ■ 

x£X xeF^ 

We deduce from (ESJ and (|SS|) : 

< sup [ip{x) + Sip{x)'j — [ sup fix) + S sup ip{x)\ < 6[ sup ip{x) — sup ip{x)\ 

xeX x£X xGFq xS-Fe x£Fo 

Since ip is Lipschitz continuous, we have 

sup ip{x) — sup ■0(2;) = sup inf (ipix) — i'iu)) 
xeFe xeFo xeF^ v^^o 

< sup inf L^\\x - y\\ 

x£F, y^Po 

= sup d{x, Fq) 

xeF, 



□ 



Corollary 20. Let X be a compact convex subset o/R". We consider an u.s.c. 
and strongly concave function ip : X W with modulus c > and a Lipschitz 
continuous function ip : X ^ M. with Lipschitz constant L.^ with respect to the 
euclidean norm. Then the maximum of (p on X is attained at a unique point 
Xq G X i.e. argmaxjf ip = {xq} and 



sup (p{x) + S^p{x)) - {p{xo) + Si;{xo)) I < L^S\ /^^^ 
xex V c 



'2(5 A/ 

[p[x) + dilj{x)) ~ [p(xo) + dil;(xo))\ < L^6\j 

x£X 

where M ~ swp^^x ''Pi^) ^ infxex ipi^). 

Proof. Define $(a;) = V'(^o) — 'p{x) for a; G X and $(a;) = +00 elsewhere. We have 
$(x) > for all X S K" and $(xo) = 0. Since 4> is l.s.c. and convex on R", then 
£ d^{xo). Moreover $ is strongly convex with modulus c. Then, using Theorem 
6.1.2 of IHUL93I Chapter VI] we have for all x,x' eX 

^x) > <^>{x') + (s I a; - x') + ^W^ - x'\\l Vs G <9$(a;') . 
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Taking x' = xo and s = we obtain for all a; G X 
which implies 

<^(x) < if{xo) - - X0II2 Vx G X . 

Using the notations of the previous Lemma, we get easily (see also Proposition 

4.32 of ESnni) for all x G F^, d{x,Fo) < where e = <5( sup^g^ ^/^(x) - 

infj;gx'0(a;))- □ 

Remark 21. To have an error estimate of the approximation of the matrix Knjt 
by the matrix Khjl, we apply Lemma |2UI in the case where 

^{x) = 'Wi{x) + Zj{x) and ^(.t) = H{x,\'wi{x)) , 

for a suitable choice of the finite elements Wi and test functions zj . Using Assump- 
tions (Hi) and {H2), we have that, for all x ^ X, \ip{x)\ < Mf\\Ww\\oo + M(, where 
llVwIloc = ||||Vu)||2||oo and Vw = (Vti;j)i<i<p. We deduce 

supV'- inf-fA < 2(Af/|lVw||oo • 

Moreover H{-,p) and H{x,-) are Lipschitz continuous with Lipschitz constants 
Lf\\p\\2 + Li and Mf respectively. Hence, "0 is Lipschitz continuous with Lipschitz 
constant 

L^ = Lf\\Vw\\oo+Li + Mf\\D'^w,\\o^ . 

5.4. Final estimation of the error of the MFEM. We now state our main 
convergence result, which holds for quadratic finite elements and Lipschitz test 
functions. 

Theorem 22. Let X be a compact convex subset of R" with non-empty interior 
and X ~ X + B2(0, where L > 0, c > 0. Choose any finite sets of discretization 
points T C M" and f cW\ Let 

Ax = ma.x{px{T),pj^{f)). 

We make assumptions (HI) and (H2), and assume that the value function at time 
t, V* , is c-semiconvex and Lipschitz continuous with constant L with respect to 
the euclidean norm, for all t > 0. Let us choose quadratic finite elements of 
Hessian c, centered at the points x of T . Let us choose, as test functions, the 
Lipschitz finite elements Zy with constant a > L, centered at the points y of T . 
For t = 0, 6, . . . ,T , let vj^ be the approximation of w* given by the max-plus finite 
element method implemented with the approximation KH,h of Kh given by (|21|l . 
Then, there exists a constant Ci > such that 

\\vl-v^U<C^{5+^) . 

When the approximation Kn^h is replaced by Kn.h, given by ()22|l . this inequality 
becomes: 

for some constant C2 > 0. 
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Proof. Let W/i and Zh denote the complete semimodules of M^ax generated by the 
families {wx)^f^f and {zy)y,^r respectively. We index the elements of T and T by 
ii, • • • ,Xp and yi, - ■ ■ ,yq respectively. Using Corollarv ll2l we have 

l^r-l'^lloo < il + j)(sup{\\P-^-{v')-v'\\^ + \\Pn.M)-v'\\o.) 

+ max \\[S^w,]h ~ S^w.Woo 

l<i<p 

To estimate the projection error ||-Pw,i(w*) — w'Hoo, we apply LemmalTTIlfor Xh = T. 
We obtain, for t £ fg, ||-P>v,i(^^*) — ^^*||oo < cdiamXAa;. Applying Lemma [T7I we 
obtain, for t E fg, \\P^^'^{v*) — f*||oo < n(a + L)Ax. Finally, using Lemma [TSl we 
get 

K-.^||oo<Ci(<5 + ^) , 

where 

Ci> {T+ l)ma.x(^cdia.mX + 7j{a + L),^{LiMf + i/A//(diamX + -) + cAlf) 

To prove the second inequality, we use Lemma [T^ together with Remark 12 II Using 
the notation of Corollary [201 and the fact that tp ^ Wi + Zj is c-strongly convex, we 
have sup 7/>-inf V' < 2(M£ + il//c(diamX + -|)) and i^, = Le + cMf + Lfc{dia.mX + 
-). We deduce that 



\{KH.hh-{KH.^)J^\ < 2{Li+cMf+Lfc{dmmX+^))^ ^ + Mf{di&mX + ^)SVS 
for i — 1, - ■ ■ ,p and j = 1, - ■ ■ ,q. Hence, there exists C2 > such that 



V 



Ax, 



lloo <C2(V^+— ) , 



5 

when S is small enough. □ 

A variant of this theorem, with a stronger assumption, was proved in |Lak03j . 

Remark 23. When T is a rectangular grid of step h > 0, meaning that T is the 
intersection of (Z/i)" with a cartesian product of bounded intervals, we have 

PxiT) < y/nh. 

Hence, when T and T are both rectangular grids of step h, we have Ax < y/nh = 
0{h) in Theorem E3 

6. Numerical results 

This section presents the results of numerical experiments with the MFEM de- 
scribed in Section |3| Wc consider optimal control problems in dimension 1 and 
2 whose value functions are known or can be computed by solving the Riccati 
equation (in the case of linear quadratic problems). 

6.1. Implementation. Wc implemented the MFEM using the max-plus toolbox 
of Scilab |Plu98| (in dimension 1) and specific programs written in C (in dimension 
2). We used the approximation KH,h of the matrix Kh- The matrix M/i can always 
be computed analytically. In all the examples below, the Hamiltonian iJ, and so 
the stiffness matrix Kh^h, have been computed analytically. We avoided storing 
the (full) matrices Mh and KH,h when the number of discretization points is large. 
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6.2. Examples in dimensionl. The next two examples are inspired by those 
proposed by M. Falcone in jBCD97j . 

Example 24. We consider the case where T = 1, = 0, X = [—1, !],[/= [0, 1], 
i{x,u) = X and f{x,u) = —xu. Assumptions {H\) and {H2) are satisfyied. The 
optimal choice is to take u* = Q whenever .t > and to move on the right with 
maximum speed {u* = 1) whenever x < 0. For all t G [0,T], the value function is: 

^ \xt if a; > 

' [a;(l--e~*) otherwise. 

We choose quadratic finite elements Wi of Hessian c centered at the points of the 
regular grid (ZAa;) n [—2,2] and Lipschitz finite elements Zj with constant a > 1 
centered at the points of the regular grid (ZAx) n X. We represent in Figure |3 the 
solution given by our algorithm in the case where 5 = 0.01, A.t = 0.005, a = 1.5 
and c = 1. We obtain a Loo-error of order 10~^. 



Figure 3. Max-plus approximation fExample 124(1 



Example 25. We consider the case where T = 1, $ = 0, X = [-1, l],U = [-1, 1], 
£{x,u) = —3(1 — |a;|) and f{x,u) = u{l — \x\). It is clear that £ and / arc bounded 
and Lipschitz continuous functions. The optimal choice is to take u* = —1 whenever 
X > and u* = 1 whenever a; < 0. Therefore, all the trajectories lie in X. For all 
t e [0, 2^] J the value function is: 

v{x,t) ==-3(1- |a;|)(l-e~*) 

We choose quadratic finite elements Wi of Hessian c and Lipschitz finite elements 
Zj with constant a. We represent in Figure 0] the solution given by our algorithm 
in the case where S = 0.02, Ax = 0.01, a ~ 2 and c = 8. Wc obtain a Loo-error of 
order 7.66 • 10"^ 

Example 26 (Linear Quadratic Problem). We consider the case where C/ = M, 

X = R, 

£{x, u) = — 7:(a:^^ + u^), f{x, u) = u, and ^ = . 
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Figure 4. Max-plus approximation (Example |23 

2 2 

The Hamiltonian is H{x,p) = ^ \ + This problem can be solved analytically. 
For a; G X, the value function at time t is 

v{x,t) = — — tanh(i)a;^. 

The domain X is unbounded and I and / are unbounded and locally Lipschitz 
continuous. We will restrict X to the set [—5; 5] so that £ and / satisfy Assumptions 
{HI) and {H2). 

We choose quadratic finite elements Wi and Zj of Hessian c = 1, centered at the 
points of the regular grid (ZAx) n [—6,6]. We represent in Figure the solution 
given by our algorithm in the interval [—1; 1] in the case where T = 5, (5 = 0.5, 
Aa; = 0.05 and L = 1. We obtain a Loo-error of 4.54 • 10~^. 




Figure 5. Max-plus approximation of a linear quadratic control 
problem (Example [SEI) 
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Example 27 (Distance problem). We consider the case where T = 1, = 0, 

X =[-!,!],[/ =[-1,1], 

iM = h -^(-I'l)' and f{.,u)^h -^(-I'l)' 
^ ^ \0 if 2;e{-l,l}, ' [0 if a;G{-l,l}. 

Putting £ = and / = on dX keeps the trajectories in the domain X but we 
loose the Lipschitz continuity of £ and /. For x £ X, the value function at time t 
of this problem is 

v{x, t) = max(— t, l^l — 1). 

Consider first quadratic finite elements Wi and Zj of Hessian c, centered at the 
points of the regular grid (ZAx) fl + Boo(0, ^)). In FigurelHl we represent the 
solution given by our algorithm in the case where 6 = 0.02, Ax — 0.01, c = 2 and 
L — 1. Since II'^'' is a projector on a subsemimodule of the Mmin-semimodule of 
c-semiconcave functions, and since the solution is not c-semiconcave for any c, the 
error of projection ||n^''(t;*) — i;*||oo does not converge to zero when Ax goes to 
zero, which explains the magnitude of the error. 




Approximated solution 
Exact solution 



.0 -0.a -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 



Figure 6. A bad choice of test functions for the distance problem fExample I27|l 

To solve this problem, it suffices to replace the test functions Zj by the Lipschitz 
finite elements with constant a > 1, centered at the points of the regular grid 
(ZAx)n[-l, 1]. This is illustrated in Figure^in the case where S = 0.02, Ax = 0.01, 
c = 2 and a = 1.1. We obtain a ioo-crror of 1.05 • 10~^. 

6.3. Examples in dimension 2. 

Example 28 (Linear Quadratic Problem in dimension 2). We consider the case 
where f/ = R^, X = R^, = 0, 

n/ \ X? + Xo U? + Uo , , 

£{x, u) = ^ — ^ — - and f[x, u) = u . 

For X £ X, the value functions at time t is 

v{x,t) = -itanh(t)(xi +xl). 
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Figure 7. A good choice of test functions for the distance problem (Example |2I|) 

As in Example the domain X is unbounded therefore £ and / do not satisfy 
Assumptions {HI) and {H2). We will restrict the domain to the set [—5; 5]^. 
We choose quadratic finite elements Wi and Zj of Hessian c centered at the points 
of the regular grid ((ZAx) n [—6, 6])^. We represent in Figurc|Hlthc solution given 
by our algorithm in the case where T = 5, 5 = 0.5, Ax = 0.1, c = 1. The Loo-error 




Figure 8. Max-plus approximation of a linear quadratic control 
problem (Example I28() 



is 9 • 10-5. 
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Example 29 (Distance problem in dimension 2). We consider the case where 

r = i,0 = o, x = [-i,i]2, (7 = [-i,i]2, 



£{x, u) 



— 1 if a; G intX, 

if X e dX, 

\ u if :r e intX, 

1 if xedX. 



For a; G X, the value function at time t is 

t) = max ( — t, max(|a;i I, \x2\) — l) • 

We choose quadratic finite elements Wi of Hessian c centered at the points of the 
regular grid ((ZAx) n [—3,3]) and Lipschitz finite elements Zj with constant a 

centered at the points of the regular grid ((ZAx) n [—1,1])^. We represent in 
Figure El the solution given by our algorithm in the case where T = 1, (5 = 0.05, 
Ax = 0.025, a = 3 and c = 1. The Loo-error is of order 0.05. 




Figure 9. Max-plus approximation of the distance problem f Example 129(1 

Example 30 (Rotating problem). We consider here the Mayer problem where T = 
l,X = 62(0, 1), U = {0}, (t){x) ^ -\xl~ Ixl, £{x, u) = and f{x, u) = {~X2, xi). 
For X Cz X, the value function at time t is 

1 3 

v{x, t) = — -(— X2sin(t) + a;iCOs(i))^ — -(a;2C0s(t) + a;isin(t))^. 

We choose quadratic finite elements Wi and Zj of Hessians Cw and Cz respectively, 
centered at the points of the regular grid ((ZAa;) n [—2,2])^. We represent in 
Figure^! the solution given by our algorithm in the case where 5 = Ax = 0.05, 
Cu, = 4 and Cz = 3. The Ltxj-error is 0.046. 
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Figure 10. Max-plus approximation of the rotating problem fExamolc 130(1 



Example 31. We consider the case where U — M., X — M.'^, (f){x) = —xf — 2x2, 

e{x,u) = ~xj — and f{x,u)^{x2,uf ■ 

We choose quadratic finite elements Wi and Zj of Hessian c^, and respectively 
centered at the points of the grids ((ZAx) n [-2,2])^ and ((ZAx) n [-11,11])^ 
respectively. We represent in Figure ITTI the solution given by our algorithm in the 
case where T = 1, 5 = 0.05, Ax — 0.025, — 10 and — 1. The Loo-error is 
0.11. (We compared the max-plus approximation with the solution of the problem 
given by the Riccati equation). 

6.4. Conclusion. We have tested our method on examples that fuUfiU the assump- 
tions of Theorem |221(see Examples 121123 EOl but also on problems that do not 
fuUfill these assumptions. The method is efficient even in the second case. The only 
difficulty comes from the full character of the matrices Mh and Kh, which limits 
the number of discretization points. To treat higher dimensional examples, we need 
higher order approximations (when the value function is regular enough). This is 
the object of a subsequent work. 
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